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:NJS  in  THE.  STAN  1  FORTH -M I TCHELL  £AI 
NUMERICAL  WEATHER  PREDICTION  CODE 


Introduction 


The  barotropic  finite  element  code  described  by  Staniforth 
and  Mitchell  in  Ref.  1  has  been  the  point  of  departure  for  exten¬ 
sive  efforts  at  the  Naval  Postgraduate  School  to  develop  improved 
techniques  for  numerical  weather  prediction.  Two  earlier  reports 
(Refs.  2  &  3)  described  applications  of  tensor  product  techniques 
to  extensions  of  the  code  developed  by  Hinsman  (Ref.  4).  The 
present  report  is  concerned  with  proposed  alterations  of  the 
Stani forth-Mi tchel 1  code  for  the  purpose  of  improving  computa¬ 
tional  efficiency. 

V. 

Three  separate  contributions  are  reported  here.  First  of 
these  is  a  scheme  for  the  direct  solution  of  the  Helmholtz  equa¬ 
tion  as  a  substitute  for  the  Concus  and  Golub  iterative  algorithm 
(Ref.  5)  used  in  the  Stani f or th-Mi tche 1 1  code.  The  second  item 
is  an  improved  solution  of  the  generalized  eigenvalue  problem 
that  arises  in  transforming  the  two-dimensional  Poisson  and  Helm¬ 
holtz  problems  into  a  series  of  one-dimensional  problems.  This 
solution  takes  advantage  of  the  special  form  of  the  matrices 
involved  (symmetric,  tridiagonal)  and  avoids  matrix  inversion. 
The  third  item  deals  with  the  special  form  of  the  generalized 
eigenvalue  problem  that  arises  when  there  is  a  periodic  boundary 
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condition  in  the  east-west  direction.  In  this  instance  the 
symmetric,  formerly  tridiagonal,  matrices  remain  symmetric,  but 


have  nonzero  elements  in  the  upper  right-hand  and  lower  left-hand 


corners.  A  separate  solution  algorithm  is  required  for  determin¬ 
ation  of  the  eigenvalues  and  also  for  the  eigenvectors. 
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The  Helmholtz  equation  may  be  written  in  the  form 


delsqw-hw=f  <1> 

where  delsq  is  the  two-dimensional  Laplace  operator,  h  is  a  func¬ 
tion  of  y,  f  is  a  function  of  x  and  y,  and  w  is  the  dependent 
variable.  For  the  Finite  Element  (FE)  discretization,  consider 
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(Note  the  periodic 
boundary  condition 
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Fig.  1.  Node  numbering  and  spacing. 


south  grid  lines,  defining  ne  intersections  (nodes).  Node  num¬ 
bering  is  from  west  to  east,  beginning  at  the  southwest  corner. 
The  positive  x  direction  is  eastward  and  the  positive  y  direction 
is  northward.  In  the  Galerkin  FE  discretization  process  the  pro¬ 
duct  h  w  is  treated  as  a  single  entity  and  the  tensor  product 
concept  is  employed  In  writing  the  matrix  form  of  the  equation. 
The  result  is 

MX  W  SY  +  SX  W  MY  -  MX  U  H  MY  =  V  <2> 

where  MX  and  SX  are  n  x  n  symmetric  "mass"  and  "stiffness"  mat¬ 
rices  for  the  x  direction,  MY  and  SY  are  the  corresponding  e  x  e 
matrices  for  the  y  direction,  U  i s  an  n  x  e  matrix  of  nodal 
values  of  the  dependent  variable,  H  is  diagonal  e  x  e,  and  V  is 
n  x  e.  MX  and  SX  depend  only  on  the  node  spacing  a(  and  MY  and 
SY  depend  only  on  the  node  spacing  b, .  Defining  equations  for 
these  matrices  are  given  in  Appendix  A.  Note  that  matrices  SX 
and  SY  are  the  negatives  of  standard  FE  stiffness  matrices.  The 


columns  of  W  contain  nodal  values  of  w  from  the  west-east  rows, 


beginning  with  the  most  southerly.  If  the  nodal  values  of  f  are 
similarly  arranged  in  an  n  x  e  matrix  F,  then  V  is  the  product 
MX  F  MY.  The  nonzero  entries  in  the  diagonal  matrix  H  are  the 
values  of  h  for  the  successive  rows  of  nodes,  beginning  with  the 
most  southerly. 

It  is  advantageous  for  the  succeeding  manipulations  to 
transpose  the  terms  of  <2>  to  obtain 

SY  WT  MX  +  MY  WT  SX  -  MY  H  WT  MX  =  VT  <3> 

where  WT  and  VT  are  the  respective  transposes  of  W  and  V.  Note 
that  all  of  the  other  matrices  are  symmetric. 

Consider  the  generalized  eigenvalue  problem 

'  SX  p,  =  1,  MX  p,  <4> 

where  p,  is  the  ith  eigenvector  (n  x  1)  and  1,  is  the  correspond¬ 
ing  eigenvalue.  Assembling  the  eigenvectors  in  an  n  x  n  modal 
matrix  P  and  the  eigenvalues  in  an  n  x  n  (diagonal)  spectral 
matrix  L,  the  result  may  be  exhibited  as 

SX  P  =  MX  P  L  <5> 

It  is  advantageous  to  require  that  the  eigenvectors  be  normalized 
so  that  PT  MX  P  =  I,  where  PT  is  the  transpose  of  P  and  1  is  the 
nth  order  identity  matrix.  If  <3>  is  pos tmu 1 t i p 1 i ed  by  P  and  <5> 
is  used  to  replace  SX  P  in  the  second  term,  the  result  is 

SY  WT  MX  P  +  MY  WT  MX  P  L  -  MY  H  WT  MX  P  =  VT  P  <6> 
Let  Q  =  WT  MX  P  and  substitute  in  <6>  to  get 

SY’  Q  +  MY  B  l  =  U  <7> 

where  SY’  =  SY  -  MY  H  and  U  =  VT  P.  Observe  that  SY’  is  not 
symmetric,  but  is  tridiagonal.  Now  the  ith  column  q,  of  Q  may  be 
found  by  solving 
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where  u,  is  the  ith  column  of  U.  The  solutions  for  the  n  q,  may 
be  assembled  into  Q.  It  is  easy  to  show  that 

U  =  P  QT  <9 > 

where  QT  is  the  transpose  of  Q.  This  completes  the  solution. 

Genera  1 i zed  Ei genva 1 ue  Problem  for  Svmmetr i c  Tr i d i aaona  1  Matr ices 

Solution  of  the  generalized  eigenvalue  problem  <4>  for  the 
symmetric  tridiagonal  ■  matrices  SX  and  MX  is  based  on  the  Sturm 
sequence  property  and  the  method  of  bisection.  The  theory  is 
given  by  Wilkinson  in  Ref.  6.  Using  this  theory  and  generalizing 
an  ALGOL  program  constructed  for  the  standard  eigenvalue  problem 
<Re.f.  7),  SUBROUTINE  TRIE1G  has  been  written  to  find  the  eigen- 
va 1 ues . 

For  any  eigenvalue,  the  corresponding  eigenvector  can  be 
found  by  rewriting  <4>  in  the  form 

(SX  -  1,  MX)  q,  •  0  <10> 

If  the  first  component  of  q,  is  chosen  as  unity,  the  first  scalar 
equation  of  <10>  determines  the  second  component  and  each  suc¬ 
ceeding  scalar  equation  determines  an  additional  component.  The 
vector  thus  found  may  then  be  normalized.  SUBROUTINE  EIGVEC  per¬ 
forms  these  operations. 

Once  eigenvectors  have  been  determined,  SUBROUTINE  RAYLEE 
may  be  used  to  refine  the  eigenvalues  by  evaluating  the  Rayleigh 
quotient  for  each  eigenvector.  Following  this  an  additional  call 
to  SUBROUTINE  EIGVEC  enhances  the  accuracy  of  the  eigenvectors. 
Experience  to  date  indicates  that  further  iterations  are  super- 


f  1  uous . 


Genera  1 i zed  Ei genva I ue  Prob 1 em  with  Periodic  Boundary  Cond i t i ons 
When  there  is  a  periodic  boundary  condition  in  the  x  direc¬ 
tion,  the  tridiagonal  form  of  matrices  SX  and  MX  is  modified  by 
the  presence  of  nonzero  elements  in  the  upper  right-hand  and 
lower  left-hand  corners.  The  matrices  remain  symmetric,  but  the 
Sturm  sequence  and  bisection  procedure  described  above  requires 
modification.  Although  the  underlying  strategy  is  unchanged, 
different  software  is  required.  In  this  instance,  the  eigen¬ 
values  are  found  by  calling  SUBROUTINE  PEREIG.  To  determine  the 
eigenvectors  SUBROUTINE  EIGVCP  is  called.  A  new  problem  arises 
here,  because,  if  the  grid  spacing  is  uniform  in  the  x  direction, 
there  will  be  double  eigenvalues.  Specifically,  if  n,  the  number 
of  subdivisions  in  the  x  direction,  is  odd,  there  will  be  a 
single  zero  eigenvalue  and  all  of  the  others  will  be  double.  If 
n  is  even,  the  eigenvalue  of  largest  absolute  value  will  also  be 
single.  Corresponding  to  double  eigenvalues,  the  eigenvectors 
are  not  unique  and  special  procedures  must  be  used  to  construct 
vector  pairs  having  the  required  orthogonality.  SUBROUTINE  EIG¬ 
VCP  tests  for  repeated  eigenvalues  and  constructs  the  appropriate 
orthogonal  pairs  of  vectors  as  required.  For  this  periodic 
boundary  condition,  eigenvalue  refinement  via  the  Rayleigh  quo¬ 
tient  is  effected  by  calling  SUBROUTINE  RAYLYP.  An  additional 
call  to  SUBROUTINE  EIGVCP  provides  better  eigenvectors.  Further 


iteration  is  believed  superfluous. 
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Remarks 

Two  additional  subroutines  are  needed  to  utilize  the  pro¬ 
posed  alterations  to  the  Stani f or th-Mi tche 1 1  program.  First  of 
these  is  SUBROUTINE  SETUP,  a  replacement  for  SUBROUTINE  EBVSET. 
The  second  is  SUBROUTINE  SLVBVP,  a  replacement  for  SUBROUTINE 
EBVP2D.  All  of  the  new  subroutines  are  listed  in  Appendix  A. 

The  routines  described  herein  have  been  tested  and  found  to 
perform  satisfactorily  on  a  12  x  12  grid.  The  parameter  RF  in 
TR1EIG  and  PEREIG  (currently  set  at  1.0  E-07)  may  need  to  be 
adjusted  for  larger  grids.  In  SUBROUTINE  EIGVCP  there  is  an 
additional  parameter,  SEP,  which  also  may  need  adjustment  for 
larger  grids.  It  is  currently  set  at  4.0  E-04. 

There  is  a  fundamental  problem  arising  from  the  use  of  single 
precision  arithmetic  in  solutions  to  the  Helmholtz  equation. 
This  results  from  the  fact  that  the  mean  geopotential  height  is 
approximately  three  orders  of  magnitude  greater  than  the  ampli¬ 
tude  of  a  representative  disturbance.  Comparisons  between  single 
and  double  precision  solutions  show  clearly  that  round-off  ad¬ 
versely  affects  the  former. 

Operation  counts  have  been  conducted  to  evaluate  the  potent¬ 
ial  saving  resulting  from  substituting  the  Helmholtz  direct  sol¬ 
ver  for  the  Concus-Golub  iterative  scheme.  For  a  13  x  13  grid 
and  assuming  only  one  cycle  of  iteration,  there  is  a  saving  of  10 
percent  for  each  time  step.  For  a  90  x  90  grid,  the  saving  is 
reduced  to  5  percent.  If  additional  iterations  are  required  for 
the  Concus-Golub  solution,  these  savings  would,  of  course,  be 
greater.  The  savings  resulting  from  the  proposed  eigenvalue  sol¬ 
utions  have  not  been  evaluated,  but  are  expected  to  be  substan- 
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APPENDIX  A. 


MASS"  and  "STIFFNESS"  MATRICES 


MX  =  i 
( n  =  4 ) 


2  ( a  i*  +a  i ) 
a  i 
0 

a  n 


a  i  0  a 

2(ai+a2)  a  2  0 

a  z  2  ( a  2  +a  3 )  33 

0  33  2(33  +a ^ ) 


Notes : 

1.  Matrices  SX  and  SY  are  negatives  of  those  called  stiff¬ 
ness  matrices  in  standard  FE  usage. 

2.  Forms  given  for  MX  and  SX  are  for  a  periodic  east-west 
boundary  condition.  For  a  Neumann  boundary  condition, 
the  terms  containing  a,  would  be  omitted. 
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FORTRAN  LISTINGS 


SUBROUTINE  SETUP (E1GVEC, E IGVAL, B IGE, B IGC, B IGA, S, HX, HY, HELM, DIR, 

1  FOURTH, NI ,NJ, CON) 

a################*########*##*###*###**######*#######*#***#*##*###****#* 

*  SUBROUTINE  SETUP  IS  DESIGNED  TO  REPLACE  EBVSET  WHEN  THE  PERIODIC  * 

#  BOUNDARY  CONDITION  IS  IMPOSED  ON  EASTERN  AND  WESTERN  BOUNDARIES  * 

a##################*#####*##*##*######*#*##*##*#***############*##*####* 
C 

C  AUTHOR:  R.  E.  NEWTON,  SUMMER  1986 

C 

C  ARGUMENTS 


c 

OUT  - 

El  GVEC 

- 

MATRIX  OF  EIGENVECTORS,  NI  X  NI,  IN  X  DIRECTION 

c 

E  IGVAL 

- 

VECTOR  OF  EIGENVALUES,  NI,  IN  NONDESCENDING  ORDER 

c 

c 

BIGE 

DIAGONAL  FACTOR  OF  COEFFICIENT  MATRIX  FOR  Y 
DIRECTION 

c 

c 

B  I  GC 

SUPERDIAGONAL  OF  UPPER  UNIT  TRIANGULAR  FACTOR  OF 
COEFFICIENT  MATRIX  FOR  Y  DIRECTION 

c 

c 

B  IGA 

SUBDIAGONAL  OF  LOWER  UNIT  TRIANGULAR  FACTOR  OF 
COEFFICIENT  MATRIX  FOR  Y  DIRECTION 

c 

IN  - 

S 

- 

SQUARE  OF  MAP  FACTOR,  N I  X  N J 

c 

HX 

- 

NODE  SPACING  IN  X  DIRECTION 

c 

HY 

- 

NODE  SPACING  IN  Y  DIRECTION 

c 

HELM 

- 

LOGICAL  SWITCH  -  TRUE  FOR  HELMHOLTZ  PROBLEM 

c 

c 

DIR 

LOGICAL  SWITCH  -  TRUE  FOR  D I R I CHLET  B.C.  ON 

POISSON  PROBLEM 

c 

c 

FOURTH 

LOGICAL  SWITCH  -  TRUE  FOR  FOURTH  ORDER  SOLUTION  OF 

POISSON  PROBLEM 

c 

NI 

-  • 

X-D I  MENS  I  ON 

c 

r* 

NJ 

- 

Y-D I  MENS  I  ON 

c 

c- 

NOTE: 

N  IS 

MAX ( N I , N J ) 

PARAMETER  <  N= 13 ) 

REAL  AS(N) , BS(N) , CS(N) , S (NI , NJ ) , HX(N I ) , HY(NJ ) , 

1  AR(N) , BR(N) , CR(N) , B IGA(N I , NJ ) , B I GC ( N I , NJ )  ,  B IGE(N I  ,  NJ )  ,  WU(  N )  . 

2  E I GVEC ( N I , N I ) , E IGVAL (N I ) , AM ( N ) , BM ( N ) , CM ( N ) , CON 
LOGICAL  HELM, DI R, FOURTH 

CF  =  2. 

I F ( FOURTH )  CF  =  5. 

C  SET  UP  "MASS"  MATRIX  FOR  X  DIRECTION  (SUBDIAGONAL  AM,  DIAGONAL  BM, 

C  SUPERDIAGONAL  CM).  PERIODIC  BOUNDARY  CONDITION  ASSUMED.  IF  NOT 

C  PERIODIC,  CALL  SETABC  INSTEAD. 

CALL  SETABX( AM, BM, CM, HX, CF, N I ) 

C  SET  UP  "STIFFNESS"  MATRIX  FOR  X  DIRECTION  (SUBDIAGONAL  AS,  DIAGONAL 
C  BS,  SUPERDIAGONAL  CS).  PERIODIC  BOUNDARY  CONDITION  ASSUMED.  IF 
C  NOT  PERIODIC,  CALL  SETD2  INSTEAD. 

CALL  SETD2N( AS, BS, CS, 1 . , HX, N I ) 

C  SOLVE  THE  GENERALIZED  EIGENVALUE  PROBLEM:  K  X  =  Z  M  X,  WHERE  K  AND 


C  M  ARE  STIFFNESS  AND  MASS  MATRICES,  RESPECTIVELY,  AND  X  IS  THE  EIGEN- 
C-  VECTOR  CORRESPONDING  TO  EIGENVALUE  2.  PERIODIC  BOUNDARY  CONDITIONS 
C  ASSUMED.  IF  NOT  PERIODIC,  CALL  TRIEIG  INSTEAD. 

CALL  PEREIGCEIGVEC, EIGVAL, AS,  BS, BM, AM, WU, NI ) 

C  SET  UP  "MASS"  MATRIX  (SUBDIAGONAL  AM,  DIAGONAL  BM,  SUPERDIAGONAL  CM) 
C  AND  "STIFFNESS"  MATRIX  (SUBDIAGONAL  AS,  DIAGONAL  BS,  SUPERDIAGONAL 
C  CS)  FOR  THE  Y  DIRECTION. 

CALL  SETABC( AM, BM, CM, HY, CF, NJ ) 

CALL  SETD2 (AS,BS,CS, 1. ,HY,NJ) 

I F( . NOT. HELM ) GO  TO  11 

C  FOR  THE  HELMHOLTZ  PROBLEM  THE  STIFFNESS  MATRIX  IS  MODIFIED  BY  SUB- 
C  TRACT  I NG  THE  CORRESPONDING  ELEMENT  OF  THE  MASS  MATRIX  DIVIDED  BY 
C  ( GZMN  TIMES  DELT»»2  TIMES  THE  RELEVANT  SQUARE  OF  THE  MAP  FACTOR). 

AS ( 1 )  =  0. 

BS ( 1 )  =  BS ( 1 )  -  BM(1)*C0N/S(1, 1) 

CS(1)  =  CS(1)  -  CM( 1 ) *CON/S (1,2) 

NJM  =  NJ  -  1 
DO  10  J  =2, NJM 

AS ( J )  =  AS ( J )  -  AM (J)*C0N/S(1,J-1 ) 

BS ( J )  =  BS ( J )  -  BM ( J ) #CON/S ( 1 , J ) 

CS(J>  =  CS(J)  -  CM (J)*CON/S(i,J+l) 

10  CONTINUE 

AS(NJ)  =  AS(NJ)  -  AM ( N J ) *CON/ S ( 1 , N J M ) 

BS(NJ)  =  BS(NJ)  -  BM ( N J ) *CON/S ( 1 , NJ ) 

CS ( NJ )  =  0. 

C  LOOP  OVER  EIGENVALUES  TO  CONSTRUCT  BIGA,  BIGC,  AND  BIGE 

11  DO  20  I  =  1,  NI 

EIG  =  E I GVAL ( I ) 

DO  12  J  =  1,  NJ 

AR ( J )  =  AS ( J )  +  E I G»AM ( J ) 

BR ( J )  =  BS ( J )  +  E I  G*BM ( J ) 

CR ( J )  =  CS(J)  ♦  E I G*CM ( J ) 

12  CONTINUE 

I F ( . NOT. HELM. AND. . NOT .DIR. AND . I . EQ . N I ) CR ( 1 )  =  0. 

C  FACTOR  COEFFICIENT  MATRIX 

CALL  SETTRI ( BR , CR , AR , AR , BR , CR , N J ) 

C  STORE  RESULTS  IN  BIGE,  BIGC  AND  BIGA.  NOTE  SIGN  CHANGES  FOR  BIGC 
C  AND  BIGA. 

DO  16  J  =  1,  NJ 

B I GE( I , J  )  =  BR ( J  ) 

BIGC(I,J>  =  -CR ( J  ) 

BIGA(I,J)  =  -AR ( J  ) 

16  CONTINUE 

20  CONTINUE 
RETURN 


cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

SUBROUTINE  PERE I G ( V , X, B , C , D, E , WU , N ) 


*  SUBROUTINE  PEREIG  USES  THE  METHOD  OF  BISECTION  TO  FIND  EIGENVALUES  * 

*  FOR  THE  GENERALIZED  EIGENVALUE  PROBLEM  INVOLVING  SYMMETRIC  (ALMOST)  * 

*  'TRIDIAGONAL  MATRICES  WITH  NONZERO  ELEMENTS  IN  'CORNERS’  (PERIODIC  * 

*  BOUNDARY  CONDITION).  IT  IS  ADAPTED  FROM  AN  ALGOL  PROGRAM  NAMED  * 

*  'BISECT'  WRITTEN  BY  BARTH,  MARTIN,  AND  WILKINSON  (NUM.  MATH.  9,  386-* 

*  393  (1967)).  'BISECT'  FINDS  EIGENVALUES  (STANDARD  EIGENVALUE  P ROB-  * 

*  LEM)  FOR  A  SYMMETRIC  TRIDIAGONAL  MATRIX.  THE  MATRICES  FOR  PEREIG  * 

*  ARE  INPUT  AS  VECTORS  -  C(N)  AS  THE  DIAGONAL  AND  B ( N )  AS  THE  SUB-  * 

*  DIAGONAL  OF  THE  "STIFFNESS"  MATRIX  AND  D ( N )  AS  THE  DIAGONAL  AND  E(N)« 

*  AS  THE  SUBDIAGONAL  OF  THE  "MASS"  MATRIX.  B(l)  AND  E(l)  ARE  THE  * 

*  "CORNER"  ELEMENTS  OF  THE  RESPECTIVE  MATRICES.  SUBROUTINE  EIGVCP  * 

*  FINDS  THE  EIGENVECTORS  AND  NORMALIZES  THEM  WITH  RESPECT  TO  THE  MASS  * 

*  MATRIX.  SUBROUTINE  RAYLYP  USES  THR  RAYLEIGH  QUOTIENT  TO  OBTAIN  IM-  * 

*  PROVED  EIGENVALUES  USING  THESE  EIGENVECTORS.  A  SECOND  CALL  TO  * 


*  TO  EIGVCP  EFFECTS  A  CORRESPONDING  IMPROVEMENT  IN  THE  EIGENVECTORS.  * 
************************************************************************ 

C 

C  AUTHOR: 

C 

C  ARGUMENTS 
C  OUT  - 

C 
C 

C  IN  - 

C 
C 
C 
C 
C 
C 

C  NOTE  - 
C 
C 

c - 

c 

INTEGER  N,Z, I , A , K , N1 

REAL  C(N),B(N),X(N), WU ( N ) , EPS  1 , EPS2 , RF , FI , XM I N , XMAX , 

1  X 1 ,  XU , XO , D ( N ) , E ( N ) , DC , DB , DD , DE , TMAX, TM I N , V ( N , N ) , 

2  Q,R,S,Q1,R1, R2, DR , QTEMP , ROLD , SW 
DATA  EPSl.RF/O. , l.E-7/ 

N1  =0 
Z  =  0 

C  CALCULATION  OF  XMAX,  XMIN 

DC=C(N) 

DB=ABS ( B ( N ) ) 

DD=D ( N ) 

DE=E ( N ) 

XMAX=(DC+DB)/ (DD+DE) 

XMIN= (DC-DB) / (DD-DE) 

C  EIGENVALUES  ASSUMED  NEGATIVE.  IF  EIGENVALUES  ARE  ALL  POSITIVE. 


R.  E.  NEWTON,  SUMMER  1986 


X 

C 

B 

D 

E  . 

WU 

N 


MATRIX  OF  EIGENVECTORS,  N  X  N,  NORMALIZED  WITH 
WITH  RESPECT  TO  "MASS"  MATRIX 
VECTOR  OF  EIGENVALUES  IN  NONDESCENDING  ORDER 
DIAGONAL  OF  "STIFFNESS"  MATRIX 

OF  "STIFFNESS"  MATRIX 
"MASS"  MATRIX 
OF  "MASS"  MATRIX 


SUBDIAGONAL 
DIAGONAL  OF 
SUBDIAGONAL 
WORK  VECTOR 
VECTOR  SIZE 


(=  N  I  ) 


MATRICES  ARE  FOR 
CONDITION.  NODE 


X-DIRECTION  WITH  PERIODIC  BOUNDARY 
SPACING  IS  HX. 


REPLACE  THE  TWO  PRECEDING  LINES  WITH  THE  FOLLOWING  TWO  LINES. 
XMAX  = (DC+DB) / (DD-DE) 

XM I N= ( DC-DB ) / ( DD  +  DE ) 

DO  2  I =N - 1 , 1 , -1 
DC=C ( I ) 

DB  =  ABS < B ( I ) ) +  ABS ( B (  1+1)  ) 

DD=D( I ) 

DE  =  E  ( I >+E( 1+1) 

TMAX= (DC+DB) / (DD+DE) 

TMIN  = (DC-DB)/ (DD-DE) 

EIGENVALUES  ASSUMED  NEGATIVE.  IF  EIGENVALUES  ARE  ALL  POSITIVE, 
REPLACE  THE  PRECEDING  TWO  LINES  WITH  THE  FOLLOWING  TWO  LINES. 
TMAX= (DC+DB) / (DD-DE) 

TM I N= (DC-DB) /(DD+DE) 

I F ( TMAX . GT . XMAX ) XMAX  =  TMAX 
I F ( TM I N . LT . XM 1 N ) XM I N=TM I N 
CONTINUE 

SET  EPS2 

I F ( XM I N+XMAX . GT . 0 . ) THEN 
EPS2=RF»XMAX 

ELSE 

EPS2=RF# ( -XM I N ) 

END  IF 

I F (EPS1 . LE. 0. )EPS1=EPS2 
EPS2=0.5*EPSl+7. *EPS2 

INNER  BLOCK 

XO=XMAX 
DO  4  I  *  1 , N 

X (  I  ) =XMAX 
WU ( I ) =XM I N 
CONTINUE 

LOOP  FOR  K-TH  EIGENVALUE 


STURM’ 


100  K=N, 1, -1 

XU=XM I N 

DO  6  I =K, 1 , -1 

I F ( XU. LT. WU ( I ) ) THEN 
XU=WU( I ) 

GO  TO  10 
END  IF 
CONTINUE 

I F ( XO. GT. X ( K ) ) XO=X ( K ) 

XI = ( XU+XO ) /2. DO 
Z  =  Z+ 1 

S  SEQUENCE 
A  =  0 

Q=C( 1 ) -X1*D( 1 ) 
Q1=C(2)-X1»D(2) 

R=B ( 1 ) -XI *E ( 1 ) 

S=C(N) -X1»D(N) 


X 


R1=B(2)-X1*E<2) 

R2=0. 

DR  =  0. 

DO  30  I *2, N-2 

IF(Q. EQ.O. ) THEN 
SW=1. 

IF (Qi+2. «R1 .EQ.O. >SU*-1. 
Q=Q1+SW*2. *R 1 
Rl =R1 +SW*Q1 

R2=SW«(B( 1+1 )-Xl#E< 1*1 > ) 
END  IF 

! F<  Q. LT. 0. )A=A*1 
S*S-R«R/Q 
ROLD=R 
R=DR-Ri«R/Q 
QTEMP*Q1-R1*R1/Q 
R 1  =  B  < I  ♦  1 )  -  X 1  *  E  ( I ♦ 1 ) -Rl »R2/Q 
Q1»C< 1*1 > -X1«D< I *1 > -R2»R2/Q 
DR*-R0LD»R2/Q 
R2=0 . DO 
Q*QT£MP 
CONTINUE 

l  F  ( Q . EQ . 0 . DO )  THEN 
SU*1 . DO 

I F ( Q1 +2. D0*R1 . EQ. 0 . DO ) SW» - 1 . DO 
Q*Q1 *2. D0*SU*R1 
Rl *R1 *SW*Q1 
R*R*SU«<B<N) -X1*ECN) > 

I  F  ( Q. LT. 0. D0)A*A*1 


ELSE 


I F(Q. LT. 0.  )  A*A* 1 
END  IF 
S*S-R»R/Q 

R»B<N)  -X1»E(N) -R*R1/Q 
Q1»Q1-R1*R1/Q 

IF<Q1.  EQ.O.  AND.  R.NE.O.  )  THEN 
A*A»  1 

ELSE 

IF(Ql.LT.O.  )  A  *  A  ♦  1 
I F  <  S-R*R/Q1 . LT. 0.  >A=A+1 
END  IF 

I F  <  A . LT . K ) THEN 

IF< A. LT. 1 ) THEN 
UU< 1 ) *X1 
XU  =  X1 


THEN 


ELSE 


UU< A+l ) *Xl 
XU*X1 

1 F ( X ( A ) . GT . X 1 ) X ( A ) =  X 1 
IF 


ELSE 


X0  =  X1 
END  IF 

IF(XO-XU. GT. 2. »RF* < ABS< XU) *ABS(X0) ) ♦EPSl >G0  TO  20 
X(K)=<X0*XU>/2. 

CONTINUE 

CALL  EIGVCPCV, B,C,D,E,X,N) 


•I  TV  TV 


-- 


CALL  RAYLYP ( X , V , B , C , WU , N ) 

CALL  E I GVCP (V,B,C,D,E,X,N) 

RETURN 

END 

cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

SUBROUT  I NE  E I GVCP (V,B,C,D,E,X,N) 


SUBROUTINE  EIGVCP  FINDS  EIGENVECTORS  BY  DIRECT  SOLUTION  OF  THE  GOV-  * 
ERNING  LINEAR  ALGEBRAIC  EQUATIONS.  FOR  ANY  PAIR  OF  EQUAL  EIGEN-  * 
VALUES,  THE  CORRESPONDING  VECTORS  ARE  MADE  ORTHOGONAL  WITH  RESPECT  * 
TO  THE  "MASS”  MATRIX.  ALL  VECTORS  ARE  NORMALIZED  WITH  RESPECT  TO  » 
THE  "MASS"  MATRIX.  * 

«**«««*«*«**«««»*««*«»«**««*«*«**««***«******#***************««««****««* 
C 

E.  NEWTON,  SUMMER  1906 


AUTHOR 


R. 


ARGUMENTS 
OUT  - 
IN  - 


V 

B 

C 

D 

E 

X 

N 


MODAL  MATRIX,  N  X  N.  (COLUMNS  ARE  EIGENVECTORS. ) 

SUBDIAGONAL  OF  STIFFNESS  MATRIX 

DIAGONAL  OF  STIFFNESS  MATRIX 

DIAGONAL  OF  MASS  MATRIX 

SUBDIAGONAL  OF  MASS  MATRIX 

VECTOR  OF  EIGENVALUES 

VECTOR  SIZE  ( =N I ) 


INTEGER  J,K,N,L,N1 
PARAMETER ( N1 *  12 ) 

REAL  V(N,N) ,P(N1) ,Q(N1) ,T(N1),X(N> ,B(N) ,C(N) 
1  H,  VN,  TN,  D I  AG,  XI , X2, VTMV, TTMV, TTMT, SUM 


D(N> , E(N) 


L  =  0 
DO  30 


10 


K=  1 ,  N 
I F(L. NE. 0)THEN 
L  =  0 

GO  TO  30 
END  IF 
XI  =X  (  K  ) 

V  ( 1 ,  K )  =  1 . DO 
V ( 2,  K ) =0. DO 
T< 1 ) =0. DO 
T(2) =1 . DO 
P<1)=B(1)-X1*E(1) 

P(2)=B(2)-X1*E(2) 

VN=-(C(1)-X1»D(1) )/P(l) 

TN*-P(2)/P(1) 

DO  10  J  =2, N-l 

D I AG=C (J)-X1*D(J) 

P(J+1)=B(J+1)-X1*E<J+1) 

V(J  +  1,K)=-(P(J) «V( J -1 , K) +DI AG*V( J , K) ) /P  <  J ♦  1  ) 
T(J+1)=-(P(J)*T(J-1)+DIAG»T(J))/P(J+1) 
CONTINUE 


C  CHECK  FOR  A  DOUBLE  ROOT 


IF(K.EQ.N)GO  TO  101 
X2  =  X  <  K  + 1 ) 

CR I T=ABS ( < X2-X1 ) / ( X2+X1 ) > 

I F ( CR IT. LT. 4. E-4 ) GO  TO  22 

C  CONSTRUCT  EIGENVECTOR  FOR  SINGLE  ROOT 

101  H=-(V(N,K)-VN)/(T(N) -TN ) 

DO  11  J  =2 , N 

V(J,K)=V<J,K)+H«T(J> 

11  CONTINUE 

C  NORMALIZE  WITH  RESPECT  TO  MASS  MATRIX 

P(1)=D<1)#V(1,K)+E(2)«V(2,K)+E(1)*V(N,K) 

DO  12  J  =2, N- 1 

P(J)=E(J)*V(J-1,K)+D(J)»V(J,K)+E(J+1)»V(J+1,K) 

12  CONTINUE 

P(N)=E(N)«V(N-1,K)+D(N)«V(N,K)+E(1)*V(1,K) 

VTMV=0 . DO 
DO  14  J  =  1 , N 

VTMV=VTMV+P  <  J ) «V ( J , K ) 

14  CONTINUE 

VTMV- 1 . CO/ SQRT ( VTMV ) 

DO  16  J  =  1 , N 

V(  J,K)=V< J,K)*VTMV 
16  CONTINUE 

GO  TO  30 

C  CONSTRUCT  EIGENVECTORS  FOR  DOUBLE  ROOT  AND  NORMALIZE 

22  L=  1 

P(l)=D(l)«V(ltK)+E(2)«V(2,K)+E(l)»V(N,K) 

Q(1)=D<1)*T<1)+E<2)*T<2)+E(1>*T<N> 

DO  23  J  =2, N-l 

P(J)=E<J)«V(J-lfK)+D(J)»V(J,K>+E(J+l)«V(J+l,K) 
Q(  J  )  =E(  J  )  *T(  J-l )  +D<  J  )  *T(  J  )  +E<  J«-l )  »T(  J  +  l ) 

23  CONTINUE 

P(N)=E(N)*V(N-1,K)+D(N)»V(N,K)+E(1)«V(1,K) 

Q(N)=E(N)*T(N-1)+D(N)»T(N)+E(1)*T(1) 

VTMV=0. DO 
TTMV=0 . DO 
TTMT=0. DO 
DO  24  J  =  1 , N 

VTMV  =  VTMV  +  P ( J ) *V ( J , K  > 

TTMV  =  TTMV+P ( J ) »T  <  J  ) 

TTMT=TTMT+Q( J ) *T( J ) 

24  CONTINUE 
H=-TTMV/VTMV 

SUM=TTMT+2. DO»H»TTMV+H*H*VTMV 
VTMV= 1 . DO/SQRT ( VTMV  > 

SUM= 1 . DO/SQRT ( SUM ) 

H=H*SUM 
DO  26  J  =  1 , N 

V(J,K+1)-V(J,K)#H  +  T(J ) »  SUM 


V(J,K)=V<J,K)»VTMV 
26  CONTINUE 

30  CONTINUE 
RETURN 
END 
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SUBROUTINE  RAYLYP(X, V, B, C, P, N) 

ft*************************************** 

*  SUBROUTINE  RAYLYP  USES  THE  RAYLEIGH  QUOTIENT  TO  FIND  IMPROVED 

*  EIGENVALUES  FROM  THE  ALREADY  NORMALIZED  EIGENVECTORS 

ft*************************************** 


C  AUTHOR 
C 

C  ARGUMENTS 
C  OUT  - 


R.  E.  NEWTON,  SUMMER  1986 


VECTOR  OF  N  EIGENVALUES  IN  NONDESCENDING  ORDER 
MODAL  MATRIX,  N  X  N.  (NORMALIZED  WITH  RESPECT  TO 
MASS  MATRIX.) 

SUBDIAGONAL  OF  STIFFNESS  MATRIX 
DIAGONAL  OF  STIFFNESS  MATRIX 
WORK  VECTOR 
VECTOR  SIZE  ( =  N I ) 


INTEGER  J , K, N 

REAL  V(N,N) ,P(N) ,X(N) ,B(N)  ,  C ( N )  ,X1 
DO  20  K= 1 , N 

P(1)=C(1)*V(1,K)+B(2)#V<2,K)+B(1)*V(N,K) 

DO  12  J  =2, N- 1 

P(J)=B(J)*V(J-1,K)+C(J)#V(J,K)+B(J+1)»V(J+1,K) 

12  CONTINUE 

P(N)=B(N)*V(N-1,K)+C(N)»V(N,K)+B(1)#V(1,K> 

XI =0. DO 
DO  16  J  =  1 , N 

X1  =  X1+P(J ) «  V  <  J  ,  K) 

16  CONTINUE 

X ( K ) =X1 

20  CONTINUE 
RETURN 
END 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

SUBROUTINE  SLVBVP(PHI , RHS, E I GVEC, B I GE, B I GC, B I GA, WK, DIR, HELM, 

1  N  I  ,  N  J  ) 


SUBROUTINE  IS  A  SUBSTITUTE  FOR  EBVP2D 


AUTHOR 


R.  E.  NEWTON,  SUMMER  1986 


ARGUMENTS 


Y 


PHI  -  SOLUTION 
RHS  -  RIGHT-HAND  SIDE 
EIGVEC  -  MATRIX  OF  EIGENVECTORS 
BIGE  -  INVERSE  OF  DIAGONAL  FACTOR  OF  COEFFICIENT  MATRIX 
BIGC  -  SUPERDIAGONAL  OF  UNIT  UPPER  TRIANGULAR  FACTOR  OF 
COEFFICIENT  MATRIX 

BIGA  -  SUBDIAGONAL  OF  UNIT  LOWER  TRIANGULAR  FACTOR  OF 
COEFFICIENT  MATRIX 
WK  -  WORK  AREA 

DIR  -  LOGICAL  SWITCH  -  TRUE  FOR  DIRICHLET  B.C. 

HELM  -  LOGICAL  SWITCH  -  TRUE  FOR  HELMHOLTZ  EQUATION 
N  I  -  X-D I  MENS  I  ON 
NJ  -  Y-D I  MENS  I  ON 


REAL  PHI <NI , NJ ) , RHS(NI , N J ) , E1GVEC(NI , N 1 > , B I GE ( N 1 , NJ ) , B l GC ( N I  ,  N J ) 
1  BIGACNI ,NJ) , WK(NJ, NI)  ,  DUM 
LOGICAL  DIR,  HELM 

C  TRANSFORM  RIGHT-HAND  SIDE 

CALL  MATPRCWK.RHS, EIGVEC, NJ,N I , NI , NJ , N I , N I , DUM, DUM, DUM, D I R, 

1  .FALSE.) 

C  PERFORM  FORWARD  REDUCTIONS 

0 

DO  5  I  =  1,  NI 

WK (1,1)  *  WK(1, I >»BIGE( I , 1) 

5  CONTINUE 

DO  15  I  3  1 ,  N I 

DO  10  J  =  2,  NJ 

WK ( J , I )  =  WK ( J , I ) #B I GE ( I , J )  +  WKCJ-l, I )*BIGA( I , J) 

10  CONTINUE 

15  CONTINUE 

C  PERFORM  BACK  SUBSTITUTIONS 

DO  25  I  =  1,  NI 

DO  20  J  *  NJ-1,  1,  -1 

WK ( J  ,  I )  =  WK ( J ,  I  )  +  WK ( J  +  1 ,  I ) *B IGC  < I , J  ) 

20  CONTINUE 

25  CONTINUE 

C  BACK  TRANSFORM  RESULTS 


DO  40  J  =  1,  NJ 

DO  35  I  =  1,  NI 
DUM  =  0. 

DO  30  K  =  1,  NI 

DUM  =  DUM  ♦  E  I  GVEC  <  I  ,  K  )  «WK  (  J  ,  K  ) 
CONTINUE 
PHI  (  I  , J  )  =  DUM 
CONTINUE 
CONTINUE 
RETURN 


%] 


SUBROUTINE  TRIEIG(V,X 

ft************************* 

t  SUBROUTINE  TRIEIG  USES  THE  METHOD  OF  BISECTION  TO  FIND  EIGENVALUES  * 

■  FOR  THE  GENERALIZED  EIGENVALUE  PROBLEM  INVOLVING  SYMMETRIC  TRIDIA-  * 
(  GONAL  MATRICES.  THE  ROUTINE  IS  ADAPTED  FROM  AN  ALGOL  PROGRAM  NAMED# 
’BISECT*  WRITTEN  BY  BARTH,  MARTIN,  AND  WILKINSON  (NUM.  MATH.  9,  386-# 
*  393  (1967)).  ’BISECT’  FINDS  EIGENVALUES  (STANDARD  EIGENVALUE  PROB-  « 

►  LEM)  FOR  A  SYMMETRIC  TRIDIAGONAL  MATRIX.  THE  MATRICES  FOR  TRIEIG  # 
i  ARE  INPUT  AS  VECTORS  -  C(N)  AS  THE  DIAGONAL  AND  B ( N )  AS  THE  SUB-  « 

>  DIAGONAL  OF  THE  "STIFFNESS"  MATRIX  AND  D ( N )  AS  THE  DIAGONAL  AND  E(N)» 
i  AS  THE  SUBDIAGONAL  ON  THE  "MASS"  MATRIX.  SUBROUTINE  EIGVEC  FINDS  » 
f  THE  EIGENVECTORS  AND  NORMALIZES  THEM  WITH  RESPECT  TO  THE  MASS  MATRIX* 


AUTHOR:  R.  E.  NEWTON,  SUMMER  1985 


ARGUMENTS 
OUT  - 


MATRIX  OF  EIGENVECTORS,  N  X  N,  NORMALIZED  WITH 
RESPECT  TO  THE  "MASS”  MATRIX 

VECTOR  OF  EIGENVALUES  IN  NONDESCENDING  ORDER 

DIAGONAL  OF  "STIFFNESS"  MATRIX 

SUBDIAGONAL  OF  "STIFFNESS"  MATRIX 

DIAGONAL  OF  "MASS"  MATRIX 

SUBDIAGONAL  OF  "MASS"  MATRIX 

WORK  VECTOR 

VECTOR  SIZE  (=  N I ) 


INTEGER  N, Ml , N, Z, I , A , K 

REAL  C(N) ,B(N) , X(N) ,WU(N) , EPS  1 , EPS2 , RF , F 1 , XMIN, XMAX, 

IQ, XI , XU, XO, DELT, D ( N ) , E ( N ), DC , DB , DD , DE, TMAX  ,TMIN,V(N,N)  ,  G ( N,  N ) 
DATA  EPS1 , RF/O. , 1 . E-7/ 

CALCULATION  OF  XMAX,  XMIN 


B( 1 ) =0. 

E( 1 ) =0. 

DC=C ( N ) 

DB=ABS ( B ( N ) ) 

DD=D ( N ) 

DE=E ( N ) 

XMAX= (DC+DB) / (DD+DE) 

XM  I  N= ( DC-DB ) / ( DD-DE ) 

NEGATIVE  EIGENVALUES  ASSUMED.  IF  EIGENVALUES  ARE  ALL  POSITIVE, 
REPLACE  THE  TWO  PRECEDING  LINES  WITH  THE  FOLLOWING  TWO  LINES. 
XMAX= (DC+DB ) / (DD-DE) 

XMIN= (DC-DB) / (DD+DE) 

DO  2  I =N- 1 , 1, -1 
DC  =C (  I  ) 

DB  =  ABS ( B ( I ) )  + ABS ( B (  1  +  1)  ) 

DD=D ( I ) 

DE=E ( I ) +  E( I +1 ) 

TMAX= (DC+DB) / (DD+DE) 


TM l N= ( DC-DB ) / ( DD-DE  > 

NEGATIVE  EIGENVALUES  ASSUMED.  IF  EIGENVALUES  ARE  ALL  POSITIVE 
REPLACE  THE  TWO  PRECEDING  LINES  WITH  THE  FOLLOWING  TWO  LINES. 
TMAX  = ( DC ♦DB )/( DD-DE ) 

TM I N= ( DC-DB ) / ( DD+DE ) 

I F (TMAX. GT. XMAX)XMAX=TMAX 
IF(TMIN.LT.XM I N ) XM I  NOTHIN 
CONTINUE 

SET  EPS2 

I F (XMIN+XMAX. GT. 0. ) THEN 
EPS2=RF*XMAX 

ELSE 

EPS2=RF* ( -XM I N ) 

END  IF 

1 F ( EPS1 . LE. 0 . ) EPS  1 =EPS2 
EPS2=0.5*EPSl+7. *EPS2 

INNER  BLOCK 

XO=XMAX 
DO  4  I =N, 1 , -1 
X( I )=XMAX 
WU( I ) =XM I N 
CONTINUE 

LOOP  FOR  K-TH  EIGENVALUE 

DO  100  K=N, 1,-1 
XU*XM I N 
DO  6  I=K, 1, -1 

I F ( XU. LT . WU  < I ) ) THEN 
XU=UU( I ) 

GO  TO  10 
END  IF 

6  CONTINUE 

I F ( XO. GT . X ( K ) )XO=X(K) 

Xl= (XU+XO) /2. 

Z=Z+1 

STURM'S  SEQUENCE 

A  =  0 
Q=l. 

DO  30  I  =  1 , N 

I F ( Q . EQ . 0 . ) THEN 

DELT=ABS ( B ( I ) -X 1 *E ( I ) ) /RF 

ELSE 

DELT= ( B ( I ) -X1»E< I > ) *»2/Q 
END  IF 

Q=C< I ) -X1«D( I ) -DELT 
IFIQ.LT.O.  ) A  =  A+  1 
)  CONTINUE 

I F ( A. LT. K ) THEN 

I F  < A. LT. 1 ) THEN 
WU( 1 )=X1 


XU  =  X1 


ELSE 

WU ( A+ 1 ) =X1 

XU  =  X1 

I F ( X ( A ) . GT .X1)X(A)=X1 
END  IF 

ELSE 

X0  =  X1 

END  IF 

IF(XO-XU.GT. 2. *RF* (ABS(XU) +ABSIXO) ) +  EPS1 )G0  TO  20 

X(K) = (XO+XU) /2. 

100  CONTINUE 

CALL  EIGVEC(V,B,C,D,E,X,N> 

CALL  RAYLEE ( V , B , C , X , N ) 

CALL  EIGVEC(V,B,C,D,E,X,N) 

RETURN 

END 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
SUBROUTINE  E I GVEC ( V , B , C , D , E , X , N ) 


*«*ft***«***tt***********ft*«****#*«ft*<t«*«*ft«**««*««***««********ft««***»»»* 

*  SUBROUTINE  El GVEC  FINDS  EIGENVECTORS  BY  DIRECT  SOLUTION  OF  THE  GOV-  * 

*  ERNING  LINEAR  ALGEBRAIC  EQUATIONS  AND  NORMALIZES  THEM  WITH  RESPECT  * 

*  TO  THE  MASS  MATRIX.  * 


AUTHOR 


R.  E.  NEWTON,  SUMMER  1985 


ARGUMENTS 
OUT  - 
IN  - 


V 

B 

C 

D 

E 

N 


MODAL  MATRIX,  N  X  N.  (COLUMNS  ARE  EIGENVECTORS.) 
SUBDIAGONAL  OF  STIFFNESS  MATRIX 
STIFFNESS  MATRIX 
MASS  MATRIX 
OF  MASS  MATRIX 
(=  N  I  ) 


DIAGONAL  OF 
DIAGONAL  OF 
SUBDIAGONAL 
VECTOR  SIZE 


INTEGER  J , K , N 

REAL  V(N,N),P(5),X(N),B(N),C(N),D(N),E(N), XI, DSQRT, SUM 
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DO  20  K= 1 , N 
XI =X ( K ) 

V  ( 1 ,  K  )  =  1 . 

P(2)=B(2)-E(2)*X1 

V(2,K)=(D(1)»X1-C(1))/P<2) 

DO  10  J  =2, N-l 

P(J+1)=B(J+1)-E(J+1)*X1 

V(J+1,K)=-(P(J)#V(J-1,K>+(C(J)-D(J)»X1)»V(J,K))/P(J+1) 

CONTINUE 


NORMALIZE  WITH  RESPECT  TO  MASS  MATRIX 


P< 1 >=D( 1 ) »V( 1, K> +E<2> »V<2,  K) 
DO  12  J  =2, N-l 


20 


12 


P(J)=E<J)»V(J-l,K)+D<J)*V(J,K)+E<J  +  l)*V(J  +  l,i<) 

CONTINUE 
P(N)=E(N)»V(N-1,K)+D(N)«V(N,K) 

SUM=0 . 

DO  14  J  =  1 ,  N 

SUM=SUM+P(J)*V(J,K) 

14  CONTINUE 

SUM  = 1 . /SQRT(SUM) 

DO  16  J  =  1 , N 

V(J,K)=V(J,K) #SUM 
16  CONTINUE 

20  CONTINUE 
RETURN 
END 

cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

SUBROUTINE  RAYLEE< V , B , C , X , N ) 


SUBROUTINE  RAYLEE  USES  THE  RAYLIEGH  QUOTIENT  TO  FIND  IMPROVED  E I  GEN- * 
VALUES  FROM  THE  ALREADY  NORMALIZED  EIGENVECTORS.  * 


AUTHOR 

R.  E.  NEWTON,  SUMMER  1985 

ARGUMENTS 

OUT  - 

X 

-  VECTOR  OF  EIGENVALUES 

IN  ASCENDING  ORDER 

IN  - 

V 

-  MODAL  MATRIX,  N  X  N. 
MASS  MATRIX) 

(NORMALIZED  WITH  RESPECT  TO 

B 

-  SUBDIAGONAL  OF  STIFFNESS  MATRIX 

C 

N 

-  DIAGONAL  OF  STIFFNESS 

-  VECTOR  SIZE  <=  N I ) 

MATRIX 

INTEGER  J , K , N 

REAL  V(N,N),P(5),X<N>,B(N),C(N>,X1 


DO  20  K= 1 , N 

P(1)=C(1)#V(1,K)+B(2)*V(2,K) 

DO  12  J  =2, N- 1 

P<J)=B<J)*V<J-1,K)+C(J)»V<J,K>+B(J+1>*V(J+1,K) 

CONTINUE 

P(N)=B(N)»V(N-1,K)+C(N)*V<N,K> 

X1=0. 

DO  16  J  =  1 , N 

X1=X1+P(J)*V(J,K) 

CONTINUE 
X  <  K  >  =X1 
CONTINUE 
RETURN 
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